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Abstract 

A strategy  to  construct  a grid  conforming  to  the  boundaries  of  a prescribed  domain  by 
using  transfinite  interpolation  methods  is  discussed.  A transfinite  interpolation  procedure 
is  combined  with  a B-spline  tensor  product  scheme  defined  by  using  suitable  control 
points.  Their  choice  is  performed  by  taking  into  account  a quality  measure  parameter 
based  on  the  condition  number  of  matrices  linked  to  the  covariant.  metric  tensors. 


1 Introduction 

The  algebraic  grid  generation  approach  relies  on  the  construction  of  a coordinate  trans- 
formation from  the  computational  domain  into  the  physical  domain.  In  particular,  this 
can  be  obtained  through  transfinite  interpolating  operators  allowing  us  the  generation 
of  grids  with  boundary  conformity.  Furthermore,  using  a Hermite-type  transfinite  in- 
terpolating scheme  we  can  obtain  orthogonal  grid  lines  emanating  from  the  boundary. 
This  can  be  very  important  for  practical  reasons  since  the  grid  point  distribution  in  the 
immediate  neighborhood  of  the  boundaries  has  a strong  influence  on  the  accuracy  of  the 
numerical  solution  of  partial  differential  equations  [5].  Furthermore,  in  case  a domain 
decomposition  is  necessary  the  orthogonality  guarantees  smoother  grids.  In  order  to  ob- 
tain a grid  with  other  specified  properties,  e.g.  the  control  of  the  shape  and  position  of 
the  coordinate  curves,  transfinite  interpolating  methods  can  be  combined  with  tensor 
product  schemes  using  suitably  chosen  control  points  (see  for  instance  [1,  2,  6,  7,  8]). 
Even  though  this  type  of  algebraic  method  is  computationally  efficient,  to  define  work- 
able meshes,  a significant  amount  of  user  interaction  is  required  for  the  selection  of  the 
control  points  involved  in  the  tensor  product.  To  overcome  this  drawback,  an  automatic 
strategy  for  choosing  the  control  points  turns  out  to  be  desirable.  Here,  following  the 
approach  first  discussed  in  [1],  we  present  an  algebraic  Hermite-type  transfinite  method 
to  construct  a grid  interpolating  the  boundary  and  its  normal  derivatives.  In  fact,  given 
a “quadrilateral”  domain  ft  C R2,  a transformation  G : R = [0, 1]  x [0, 1]  — > ft  is  defined 
as 

G(s,  t)  :=  TP(s , t)  + (P:  ® P2)([d>,  4>]  - TP)(s,  t)  (1.1) 

where  Tp  is  a tensor  product  surface  i.e.  Tp(s,t)  :=  Y1T=\  J2j=i  with 

Bit 3 denoting  the  usual  cubic  B-spline,  4>  and  ip  are  boundary  curves  and  (Pi  ®P2)  is  the 


2 
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Boolean  sum  of  Hermite-type  blending  function  linear  operators.  The  set  Q = {Q,j,  i = 

1 j = 1, . . . , n}  is  the  set  of  control  points. 

As  already  noted,  the  choice  of  the  control  points  is  a crucial  matter.  In  this  paper 
we  take  into  account  a grid  quality  measure  parameter  for  their  selection.  In  particular, 
the  proposed  automatic  procedure  relies  on  the  fact  that  some  grid  properties  can  be 
described  in  terms  of  the  condition  number  of  matrices  linked  to  the  covariant  metric 
tensors  [4] . Therefore,  the  control  points  are  chosen  minimizing  their  condition  number. 

The  outline  of  this  paper  is  as  follows.  In  Section  2,  the  transformation  (1.1)  is  given 
in  detail  and  its  properties  are  investigated.  In  Section  3,  a way  for  choosing  the  control 
points  is  proposed  relying  on  a particular  quality  measure  parameter.  Finally,  in  Section 
4 some  numerical  results  are  presented  to  illustrate  the  features  of  the  proposed  strategy. 

2 The  transformation 

In  this  section  the  transformation  (1.1)  is  characterized.  Let  us  consider  a “quadrilateral” 
domain  fl  c 1R2  such  that  dfl  = uf=1dfli,  with  dfli,dfl2,  dfl3,  df^  being  the  supports  of 
four  regular  curves  % : [0, 1]  — > dfl, , i = 1, . . . , 4 taken  counterclockwise.  Furthermore, 
let  us  suppose  that  dfli  fl  dfl3  = 0 and  dfl2  D dfl 4 = 0,  with  any  other  intersection 
occuring  only  at  the  end  points  of  the  boundary  curves.  In  particular,  the  following 
compatibility  conditions  are  assumed 

7i(0)  = 74(1),  7i(l)  = 72(0),  72(1)  = 73(0),  74(0)  = 73(1)  • 

For  later  convenience,  we  set  </>i(s)  :=  7i(s),  02 (s)  :=  73(1  - s ) denoting  by  s the 
curve  parameter  running  on  [0, 1]  and  we  set  0i (t)  :=  74(1  - t),  02(t)  :=  72 (t)  denoting 
by  t the  curve  parameter  running  on  [0, 1],  In  addition,  the  components  of  the  0-curves 
and  0-curves  are  denoted  by  0x,0y  and  ipx,ipy  respectively. 

Next,  we  define  four  additional  curves  by  computing  the  derivatives  of  the  <p  and 
0-curves,  i.e., 


Ms)  = ii^HtfOO)',  (4>Us)Y)  , i = 1,2, 

(2.1) 

Vb+ 2(t)  = jj^jj ;(-(V^(t))',  (ip jit))')  , 3 = 1,2  , 

with  C a constant  value  also  depending  on  the  curve  orientations  and  with  ||  • ||2  the 
Euclidean  norm.  Then,  we  introduce  the  linear  operators 

Pi[0](s,f)  :=  TA=iai(t)Ms)  > !>2[0](s.t):=Ej=iajW}(t). 

(2.2) 

PiP2[0,0](s,t)  :=  Ei=i  (ai(t)P2[ip](s,Ui)  + ai+2(t)—2^’8](SlUi)')  , 
where  = 0,  u2  = 1.  The  functions  a*,  % = 1,...,4,  are  the  dilated  versions  of  the 
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classical  Hermite  bases  with  support  on  [0,  n]  and  on  [1  - u,  1]  being  0 < u < 1,  i.e. 

:=  (1  + 2t)(l  - f)2  > q3(s)  :=  s(l  - f )2  , s€[0,«], 

a2(S):=(3-2^1)(^l)2  , a4(s)  :=  (s  - 1)(^)2  , se[l~u,l}. 

(2.3) 

The  Boolean  sum  operator  (Pi©P2)  — Pi  + P2  — P1P2  provides  the  blending  function 
surface 

B(s,t)  :=  (P1  ® P2)[0,0]M)  = P1ms,t)  + P2[0](s,t)  - PiP2[4>M(s,t)  . (2.4) 

It  is  known  that  B satisfies 

B(uit  t)  = Ut)  , < = 1, 2 , i = 3, 4 , . , 

B{s,Wj)  -<j>j(s)  , j = 1,2  ° ar1-  = <t>j(s)  , j = 3,4  , 

where  ui  = U3  = 0,  u2  = U4  = 1 and  rci  = 1^3  — 0,  w2  = 104  = 1.  It  is  worthwhile  to 

remark  that,  as  we  are  dealing  with  orthogonal  grid  lines  emanating  from  the  boundary 
of  the  domain,  the  intersecting  boundary  curves  must  be  also  orthogonal.  Thus,  the 
following  additional  conditions  are  assumed: 

<f>i+ 2(0)  = <j)i+2{  1)  = lp'2{Wi)  , 

^+2(0)  =4>\ (Ui),  ipi+ 2(1)  =4> 2(u,-)  , * = 1)2.  (2.6) 

4>'i  (0)  = V’i  {wi ) , 4>"  (1)  = i>2  (u>i ) , 

Now,  in  order  to  define  a suitable  grid,  following  the  approach  given  in  [1],  we  use 
the  linear  transformation  G 

G(s,  t ) :=  TP(s,  t ) + (A  ® P2)  ([0, 0]  - TP)  (s, t)  (2.7) 

where  Tp(s,t)  :=  I2j=i  QijBlfi(s)Bj^(t)  with  Pl  3 denoting  the  usual  cubic  B- 

splines  with  uniform  knots.  The  set  Q = {Q,j,  i = 1, . . . ,m,  j = I, . . . , n}  is  a suitable 
set  of  control  points  whose  definition  is  discussed  in  Section  3.  It  should  be  noted  that 
in  (2.7)  the  Boolean  sum  operator  is  also  acting  on  a surface  Tp(s,t).  In  this  case 
(2.2)  is  used  taking  the  eight  boundary  curves  Tp(0,  t),  Tp(l.t),  Tp(s, 0),  Tp(s,  1), 

dTp(0,t)  dTp(l,t)  dTp(sfi)  dTp(s,  1) 

ds  ’ ds  * dt  * Ot 

It  is  easy  to  show  that  G still  satisfies  G(ui,t)  = ipi(t)  , i = 1,2,  , i = 

3,4,  G(s,Wj)  = <f>j(s)  , j = 1,2  and  dG^~i)  = cpj(s)  , j = 3,4.  Furthermore,  because 
of  the  locality  of  the  blending  functions  a,,  i = 1, ....  4,  the  control  of  the  coordinate 
lines  obtained  by  means  of  the  evaluation  of  G over  a parameter  set  in  the  interior 
of  the  domain  is  mainly  based  on  the  contribution  of  Tp.  This  fact  and  the  use  of  B- 
splines  ensures  the  convex-hull  property  in  the  interior  of  the  domain.  This  property  is  of 
importance  in  numerical  grid  generation  to  locate  the  grid  with  respect  to  the  position 
of  control  points. 
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3 Grid  quality  measure 

It  is  well  known  that  grid  generation  techniques  sensible  to  grid  quality  features  are 
particularly  attractive.  Thus,  in  this  section,  we  discuss  a strategy  to  choose  the  set  Q 
of  control  points  based  on  a suitable  grid  quality  measure  parameter. 

Given  a set  of  grid  points  Q :=  defining  the  quadrilateral  cells  {Ci 

quality  measures  can  commonly  include:  grid  “skewness”,  measuring  the  departure  of 
Cij  from  a rectangle,  grid  “aspect  ratio” , measuring  the  departure  of  Cij  from  a rhombus 
or  grid  “conformality”,  measuring  the  departure  of  C \j  from  a square  (see  for  instance 

[5])- 

Here,  as  done  in  [4]  for  the  case  of  unstructured  grids,  we  define  a grid  quality  measure 
taking  into  account  the  condition  number  of  particular  matrices  derived  from  the  grid. 
As  explained  below,  somehow  this  quality  parameter  measures  the  departure  of  Cij  from 
a square. 

The  strategy  starts  with  a set  Q of  control  points  obtained  by  evaluating  on  a coarse 
parameter  set  Sc  = a Lagrange  blending  function  surface  (for  detail 

related  to  Lagrange  blending  function  methods  we  refer,  for  instance,  to  [3])  by  working 
only  with  the  four  boundary  curves  of  the  given  domain.  Then,  using  Q1  a first  grid  is 
obtained  by  evaluating  the  surface  G in  (2.7)  on  a fine  parameter  set  Sf  = {(sj,fy)}^’^ l 
obtaining  the  grid  points 


g :=  {Gij  = = G(si,tj),  i = 1,...  ,M,  j = 1,.,.  ,N) 


The  set  Q is  then  used  to  define  (M  - 1)  x (N  - 1)  bidimensional  matrices  associated 
with  the  ( M — 1)  x (N  - 1)  quadrilateral  cells  , i = 1, . . . , M — 1,  j = 1, . . . , N — T. 
These  matrices  are  defined  as 


Ai 


hi 


Ui+ l,j 

GSVij 


Gx 

i,j 

- Gli 


rtx  rtx 

Gri,i+ 1 
r*y  _ nv 


) , i = 1,  • • • , M — 1,  j = 1, ...  ,N  — 1 


(3.1) 


and  their  condition  number  K(Aij)  is  related  to  the  stretch  of  the  cells.  In  fact,  it  is  easy 
to  prove  that  K(Aij)  :=  | 1 12  • \ \A~Jj  H2  = 1 if  and  only  if  we  are  dealing  with  a cell 
Cij  where  the  three  points  Gij+\ , Gtj,  generate  half  a square  [9].  On  the  other 

hand,  in  order  to  involve  all  the  grid  points  in  the  quality  measure  it  is  also  convenient 
to  define  the  boundary  matrices 


1 


A 


M-l,j 


/~1X 

LTi+l,JV 

(*i+l,N 


fix 

UM,j+l 

r<v 

"Afj+l 


fix  fix  fix  \ 

i,N  'Jri,N—l  ^ i,N  \ 

_ fiy  fiy  _ fiy  h 
^ i,N  ^i,N- 1 i,N  / 

fix  fix  fx  \ 

~UM,j  \ 

_ f<y  fiy  —ny  I 

M,j  ^ M,j  / 


i = 1,. . . ,M  - 1, 

, j = 1,...,N  -1, 


(3.2) 


A 


Gx  ctx  rtx  /~ix  \ 

M,N  " ^ M,N  ^M,N- 1 \ 

rty  _ r*y  r<y  _ r*y  I 

^M,N  ^ M,N  J 


so  that  the  boundary  points  are  also  taken  into  account. 

Next,  we  modify  the  initial  set  Ql  of  control  points  minimizing  the  following  objective 
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function 


u = ih  (e^t1  eJL_ix  K(Aitj) + Eflr1  k{a\jN_  x)+ 
EjLl1  K{Alt_hj)  + KiAl'^-S)  • 


(3.3) 


The  minimization  is  done  with  respect  to  the  control  points  under  suitable  constraints 
on  their  coordinates  depending  on  the  geometry  of  the  domain  Cl.  This  is  the  only  user 
interaction  required. 

Obviously,  since  ideal  inner  cells  are  characterized  by  an  associated  matrix  Aij  having 
a condition  number  close  to  one,  the  optimal  distribution  of  the  control  points  should 
guarantee  ming  f0b  ~ 1.  On  the  other  hand,  ming  fDb  strongly  depends  on  the  geometry 
of  the  domain  (for  example  in  case  of  a squared  domain  the  optimal  value  is  ming  f0b  = 1 
while,  in  general,  this  value  is  not  reached). 


Summary  of  the  Method 

(1)  Compute  the  initial  set  of  control  points  Q1  by  means  of  a Lagrange  blending  function 
method  using  the  four  given  boundary  curves, 

(2)  Compute  the  initial  grid  Q%  = {G(sj,  tj),  i = 1, . . . , M,  j = 1, . . . , N}  with  G given 
in  (2.7)  by  using  the  set  of  control  points  Q', 

(3)  Minimize  the  objective  function  (3.3)  so  defining  a new  set  of  control  points  Q? , 

(4)  Compute  the  final  grid  Qs  = {G(Si,tj),  i — 1, . . . ,M,  j = 1, . . . ,N}  with  G given 
in  (2.7)  by  using  the  set  of  control  points  Q s with  M » M,  N » N. 

Remark  3.1  We  note  that,  in  order  to  reduce  the  computational  cost  of  the  minimiz- 
ation procedure,  the  integers  M and  N are  chosen  less  than  M and  N. 


4 Numerical  Results 

We  conclude  the  paper  giving  some  numerical  results  testing  the  properties  of  the  trans- 
formation G and  showing  the  performance  of  the  proposed  approach. 

Three  domains  are  considered.  For  each  of  them  we  present  the  initial  grid  obtained 
by  the  transformation  G using  the  initial  set  of  control  points  Q 1 and  the  final  grid 
obtained  using  the  set  of  control  points  Qf  resulting  from  the  minimization  procedure. 
In  all  the  figures  the  control  points  are  denoted  by  the  symbol  **’.  The  minimization 
problem  is  solved  by  using  a sequential  quadratic  programming  method  i.e.  by  using  the 
routine  constr  of  the  Optimization  toolbox  of  the  Matlab  package.  In  the  minimization 
procedure,  the  constraints  on  the  control  points  are  chosen  so  that  some  geometric 
properties  of  the  domain,  such  as  symmetry  and  convexity,  are  preserved.  Furthermore, 
in  all  the  examples  M and  N are  equal  to  M and  to  N.  The  values  of  the  objective 
function  before  the  minimization  (flb)  and  after  the  minimization  (//h)  are  also  given  in 
the  figure  captions. 

The  first  and  the  second  test  display  a “waterway”  grid  and  a "W-shaped  grid  with 
their  control  points  before  and  after  the  minimization  procedure.  The  effectiveness  of 
the  method  is  evident. 
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Figure  3 shows  a grid  composed  of  six  sub-grids,  obtained  via  a domain  decompos- 
ition approach.  In  this  case,  the  Hermite-type  interpolation  method  guarantees  a C1 
connection  among  the  patches.  Here,  the  two  values  of  foh  and  fjth  in  the  figure  captions 
refer  to  the  “horizontal”  and  “slanted”  grids,  respectively. 
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